DGE post SCI and MSC Treatment
Ramil Hakim
07 augusti 2017
1 Project Description
Hypothesis: Transplantation of syngeneic bone-marrow derived mesenchymal stem cells (BMSCs) has superior therapeutic effect in comparison to transplantation of allogenic BMSCs following experimental spinal cord inury (SCI) in mouse.
Experimental setup & sequencing: C57BL/6J female mice underwent a severe (75 kdyn) contusion SCI followed by transplantation of bMSCs 24 h post SCI. Animals were sacrificed 20 days following SCI. Total RNA, digested of DNase, was isolated and sequenced (125 cycles paired-end) in one lane using the HiSeq2500 system and v4 sequencing chemistry (Illumina Inc.) performed by the SNP&SEQ Technology Platform (Uppsala, Sweden).
Data analysis: A read count matrix was obtained from the sequencing core facility. Data was analyed using the edgeR and limma packages (both available through bioconductor.org) using R version 3.3.2. Two additional key packages used were ggplot2 and data.table.
2 Data Overview
Table 1. Dataset characteristics
| Characteristic | Value |
|---|---|
| Samples (n): | 12 |
| Groups (n): | 4 |
| Unique ENSEMBL IDs (n): | 46078 |
3 Data Pre-Processing
3.1 Transforming from raw scale and removing lowly expressed genes
Background: The expression of a gene must reach a certain threshold for it to be translated into a protein. Translation into a protein is a prerequiste for a gene to have any biological function. Thus, genes with low number of read counts across samples are probably not differentially expressed and should be removed. However, a greater sequencing depth (i.e. a larger library size) will result in a higher read count, thus introducing a bias to the DGE analysis. Therefore, prior to filtering, raw counts are transformed into counts per million (CPM) which accounts for the difference in library size.
Transformation: Raw counts are transformed into CPM by dividing each count with the sum of the column (i.e. the library size) and then multiplying by 1e6. Log-CPM is calculated by taking log2 (raw count + 0.25).
Filtering: Lowly expressed genes are removed from the count matrix by filtering. A gene is defined as highly expressed if CPM>1 for at least three samples. Three samples was chosen since it is equal to the smallest group size (which is equal between groups in this case).
Fig 1. Density of log-CPM values pre -and post filtering
Fig 1. Figure reports the density of log-CPM for every sample (by color) pre -and post filtering of lowly expressed genes. Filtering is conducted using log-CPM values. Vertical dashed line represents the cut off (log-CPM=0). The figure shows a distinct shift of the density from below the threshold (Fig 1A) to above the threshold (Fig 1B).
3.2 Normalizing gene expression distributions
Background: The read count is affected by: 1) the gene expression and 2) the sequencing depth. The sequencing depth equals the library size. The library size is defined as the sum of counts for each sample (i.e. column sum). The counts per sample represent the relative abundance of each gene. Highly expressed genes can consume a substantial proportion of the library size thus making the other genes seem underexpressed. Therefore, normalization is conducted in order to ensure that the distribution of the expression is similar for each sample. All samples should have a smiliar range and distribution of expression (log-CPM).
Normalization: Scaling factors are calculated using the trimmed mean of M-values (TMM). The algorithm finds a set of scaling factors which minimizes the log-fold change between the samples. Scaling factors >1 downscale the counts while scaling factors <1 scale the counts upwards. The effective library size is then obtained by taking the product of the original library size and the scaling factor for each sample respectively.
Fig 2. Distribution of log-CPM pre -and post normalization
Fig 2. Figure reports the distribution of gene expression (log-CPM) for each sample. Fig 2A reports the distribution prior to normalization while Fig 2B reports the distribution following normalization of library sizes using trimmed means of M-values. Boxplots are based on all log-CPM values while points represent a random sample of 1e4 observations (due to processing time issues). The difference in the distribution of log-CPM using original and effective library sizes is minor but adjusted for.
4 Unsupervised Clustering 1: PCA
Principal component analysis (PCA): A PCA aims at producing a low-dimensional representation of the dataset. The principal components are each normalized linear combinations of a set of features constructed with loadings with the intention to achieve maximal variance. Normalized means that the sum of squared loadings equals one. Furthermore, the principal components are constructed to be uncorrelated to each other.
Fig 3. Variance explained by principal components based on the 500 genes with highest variance
Fig 3. Figure reports the proportion variance explained by each principal component. Fig 3A reports the proportional variance explained by each component while Fig 3B reports the cumulative variance explained by the components. It is obvious that the first principal component explains the majority of the variance while the remaining components explain only a small portion of the variance.
Fig 4. Estimation of distribution of the variance explained by the principal components
Fig 4. Figure reports the distribution of bootstrapped proportion of variance explained by each principal component. Bootstrap is conducted using the 500 genes with the highest variance. Fig 4A reports the distribution for the first principal component while Fig 4B reports the distribution for the remaining eleven components (due to differences in magnitude). Shaded region in Fig 4A represents a bootstrapped nonparametric 95 % confidence interval. The histograms show narrow distributions which confirms the observations and conclusion in Fig 3.
Table 2. Upper and lower bounds (95 % confidence interval) for the proportion of variance explained by principal component 1 to 11
| PC1 | PC2 | PC3 | PC4 | PC5 | PC6 | PC7 | PC8 | PC9 | PC10 | PC11 | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Upper bound : | 0.775 | 0.044 | 0.038 | 0.027 | 0.022 | 0.011 | 0.009 | 0.006 | 0.004 | 0.005 | 0.003 |
| Lower bound: | 0.823 | 0.060 | 0.051 | 0.042 | 0.033 | 0.018 | 0.014 | 0.010 | 0.007 | 0.006 | 0.004 |
Fig 5. Multidimensional scaling plot of the 500 genes with the highest variance
Fig 5. Figure reports the distribution of the samples along the first and second principal component. Study groups are indicated with ellipses. Ellipses represent a 95 % confidence interval based on the samples within each group. There is a distinct separation between the four study groups. Uninjured samples separate from the three other study groups which contain animals exposed to SCI. The syngeneic and allogenic groups do form an “venn diagram like” pattern in relation to the injury only group in the sense that they both have one animal which overlaps with the injury control study group.
Fig 6. K means clustering (4 groups) of 1000 boostrap replicates for the first -and second principal component for the 500 genes with the highest variance
Fig 6. Figure reports the first and second principal component for 1000 bootstrapped replicates per sample based on the 500 genes with highest variance. K means clustering with 4 groups and 20 starts was conducted on both the bootstrapped data. Ellipses are centered around the mean value of the botstrapped values for each cluster and represent a 95 % confidence interval for the data in each cluster. The clustering confirms the hypothesis postulated in Fig 5.
Fig 7. Top 10 positive and negative loadings for the first -and second principal component
Fig 7. Figure reports the top 10 positive and negative loadings for the first -and second principal component.
Fig 8. Hierarchical clustering of samples based on 1000 genes with highest variance
Fig 8. Figure reports a dendrogram which illustrates hierarchical clustering of samples based on the 1000 genes with the highest variance. Uninjured animals and animals which received allogenic transplantation cluster together and separately. The injured only animals and one of the syngeneic animals cluster together which supports the intercept like pattern observed in the PCA.
5 Negative Binominal Dispersions by Weighted Likelihood Empirical Bayes
Fig 9. Cox-Reid profile-adjusted likelihood estimated tagwise, common and trended dispersions
Fig 9. Figure reports tagwise, common and trended dispersions.
6 Linear Modelling and Empirical Bayes Moderation Using Precision Weights
Background: The variance of raw counts is not independent of its mean. Furthermore, the variance of raw counts increases with the count size (opposite is true for log-counts).
Voom transformation: Acronym for mean-variance modelling at the observational level. The mean-variance relationship is estimated in the data which is then used for computing precision weights for each gene. The precision weights are then implemented in the linear modelling in order to adjust for heteroscedasticity. Estimation of precision weights is done as follows: 1) log-CPM=log2(raw count + 0.5), 2) a linear modelled is fitted to each gene, 3) lowess function is fitted to the scatterplot between average log-CPM and sqrt(st.dev), 4) the variance of the log-CPM for each gene is predicted using the lowess trend line 5) the inverse of the variance for each log-CPM value is the precision weight.
Empirical Bayes moderation of standard error: Ranks genes in order of evidence of differential expression.
Fig 10. Mean-variance relationship pre -and post voom transformation
Fig 10. Figure reports the mean-variance relationship pre -and post application of the voom function. Fig 9A reports the average log-CPM against the quarter root of the variance. Fig 9B reports average log-CPM against the log2(st.dev). Blue line reports the average log2(st.dev). The red line is a linear trend fitted to the black dots. Each black dot represents a gene. Fig 9A illustrates that the variance is decresing when the average expression is increasing. In Fig 9B the dependency is removed and the mean variance is unchanged when the average expression increases.
7 Differential Gene Expression
Determining significant gene differential expression for a contrast: A simple Bayesian model is used for moderating standard errors across genes (squeeze them towards a common value) and producing moderated t-statistics. These moderated t-statistics are used for significance analysis. Moderated t-statistics have a higher degree of freedom in comparison to usual t-statistics due to the increase in reliability resulting from the smoothening of standard errors. P-values are adjusted for multiple testing using Benjamini and Hochberg’s method to control for false discovery rate (FDR).
Fig 11. Number of differentially expressed genes for each contrast
Fig 11. Figure reports venn diagrams containing the number of differentially expressed genes for each contrast. Fig 10A reports the number of differentially expressed genes for study groups which received a SCI while Fig 10B reports the study groups in relation to the uninjured samples.
Table 3. Number of differentially over -and under expressed genes for each contrast
| Injury-Syngeneic | Injury-Allogenic | Injury-Uninjured | Syngeneic-Allogenic | Syngeneic-Uninjured | Allogenic-Uninjured | |
|---|---|---|---|---|---|---|
| Downregulated: | 2 | 16 | 1503 | 60 | 101 | 1314 |
| No change | 15444 | 15434 | 11206 | 15390 | 14443 | 11505 |
| Upregulated: | 6 | 2 | 2743 | 2 | 908 | 2633 |
Fig 12. Mean difference -and volcano plot
Fig 12. Figure 11A reports a mean-difference plot which illustrates the number of over -and under expressed genes. Threshold is set at log2(fold change) +/-1 (blue lines). Blue dots represents genes above or below the log-fold change thresholds while red dots represent those genes which are above/below th thresholds and are significantly differentially expressed. Fig 11B is a volcano plot which reports the number of significantly over -and underexpressed genes (marked with red). Data in both plots is for the contrast “syngeneic IDmBMSC vs allogenic IDmBMSC”.
Table 4. 10 most significantly up -and downregulated differentially expressed genes
| Gene | log2(fold change) | P-value (adjusted) | Gene | log2(fold change) | P-value (adjusted) |
|---|---|---|---|---|---|
| Gimap4 | -3.12 | 0.0014 | ENSMUSG00000102729 | 1.31 | 0.0406 |
| Cd8a | -5.07 | 0.0007 | A230065H16Rik | 1.57 | 0.0409 |
| Ms4a4b | -4.43 | 0.0025 | NA | NA | NA |
| Lck | -2.92 | 0.0029 | NA | NA | NA |
| Tbc1d10c | -3.09 | 0.0044 | NA | NA | NA |
| Cd8b1 | -6.42 | 0.0029 | NA | NA | NA |
| Tigit | -5.78 | 0.0006 | NA | NA | NA |
| Acap1 | -3.19 | 0.0044 | NA | NA | NA |
| Ctsw | -3.37 | 0.0054 | NA | NA | NA |
| Cd2 | -4.48 | 0.0044 | NA | NA | NA |
8 Gene Ontology and KEGG Enrichment Analysis
Gene Ontology (GO): A major bioinformatics initiative which offers a computational representation of the biological function of genes at molecular, cellular and tissue level. This tool enables one to annotate genes with their function.
KEGG (Kyoto Encyclopedia of Genes and Genomes): KEGG pathway are manually drawn pathway maps which represent the current knowledge within metabolism, cellular processes and many more.
Fig 13. Top 10 enrichment terms/pathways
Fig 13. Figure reports gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis. Top 10 enrichment IDs (term and pathway respectively) are reported. Horizontal axis represents the number of differentially expressed genes associated with each ID. Values indicate the p-value for each ID respectively.
9 Unsupervised Clustering 2
9.1 Agglomerative hierarchical clustering with heat map
Background: Purpose of hierarchical clustering in combination with heatmap is to find subset of genes which explain the difference between study groups.
Hierarchical clustering: An unsupervised statistical learning method which aims at clustering the data in a not pre-determined amount of clusters. In agglomerative hierarchical clustering the tree is built from the terminal nodes (leaves) towards the root. A dendrogram is utilized to displayed the clusterings. A heatmap is added to the hierarchical clustering to enable the identification of gene clusters which might explain the hierarchical clustering.
Fig 14. Hierarchical clustering of samples together with heatmap of significantly differentially expressed genes
Fig 14. Figure reports a heat map with hierarchical clustering (indicated with dendrograms) using log-CPM values. Only significantly differentially expressed genes are included (and genes with NA symbols were removed). A cluster in the top 10 rows is visible. In this cluster the allogenic IDmBMSCs cause the genes to be up-regulated, injury causes them to be somewhat up-regulated while syngeneic IDmBMSCs causes the lowest amount of up-regulation with a tendency towards down regulation. Thus syngeneic IDmBMSC is most similar to uninjured animals.
9.2 Multiscale bootstrap resampling of agglomerative hierarchical clustering
Fig 15. Multiscale bootstrap resampling of hierarchical clustering
Fig 15. Figure reports hierarchical clustering using multiscale boostrap resampling with 1000 replicates. Numbers with red color are the approximately unbiased (AU) p-values while blue numbers represent bootstrap probability (BP). AU >95 and BP >70 are considered to be significant clusters. Red boxes indicate clusters with AU >95.
Fig 16. Unsupervised clustering using the 1000 genes with the highest F value (between all contrasts) 
Fig 16. Figure reports a hierarchical clustering based on the 1000 genes with the highest F statistica (higher F statistica indicate a larger difference between contrasts). The conclusion is as for figure 8, i.e. that the injured animals do show similarites and that the ‘intercept’ between the syngeneic treated animals and the injury control is seems to ber present.
9.3 Description of function of ‘cluster genes’
Apobec3: Selectively targets single-stranded DNA and does not deaminate double-stranded DNA or single-or double-stranded RNA. Exhibits antiviral activity against HIV-1, simian immunodeficiency viruses (SIVs), mouse mammary tumor virus (MMTV) and friend murine leukemia virus (FrMLV) and may inhibit the mobility of LTR retrotransposons. [ref: UniProt]
Epsti1:
Akna: Transcription factor that specifically activates the expression of the CD40 receptor and its ligand CD40L/CD154, two cell surface molecules on lymphocytes that are critical for antigen-dependent-B-cell development. Binds to A/T-rich promoters. [ref: UniProt]
Itgal: Involved in a variety of immune phenomena including leukocyte-endothelial cell interaction, cytotoxic T-cell mediated killing, and antibody dependent killing by granulocytes and monocytes. Contributes to natural killer cell cytotoxicity. Involved in leukocyte adhesion and transmigration of leukocytes including T-cells and neutrophils. [ref: UniProt]
Rac2: Plasma membrane-associated small GTPase which cycles between an active GTP-bound and inactive GDP-bound state. In active state binds to a variety of effector proteins to regulate cellular responses, such as secretory processes, phagocytose of apoptotic cells and epithelial cell polarization. Augments the production of reactive oxygen species (ROS) by NADPH oxidase. [ref: UniProt]
H2-Q7: Involved in the presentation of foreign antigens to the immune system. [ref: UniProt]
Sash3: May function as a signaling adapter protein in lymphocytes. [ref: UniProt]
Psmb8: The proteasome is a multicatalytic proteinase complex which is characterized by its ability to cleave peptides with Arg, Phe, Tyr, Leu, and Glu adjacent to the leaving group at neutral or slightly basic pH. The proteasome has an ATP-dependent proteolytic activity. This subunit is involved in antigen processing to generate class I binding peptides. May be involved in the inflammatory response pathway. Required for adipocyte differentiation. [ref: UniProt]
Irf1: Transcriptional regulator which displays a remarkable functional diversity in the regulation of cellular responses. These include the regulation of IFN and IFN-inducible genes, host response to viral and bacterial infections, regulation of many genes expressed during hematopoiesis, inflammation, immune responses and cell proliferation and differentiation, regulation of the cell cycle and induction of growth arrest and programmed cell death following DNA damage. Stimulates both innate and acquired immune responses through the activation of specific target genes and can act as a transcriptional activator and repressor regulating target genes by binding to an interferon-stimulated response element (ISRE) in their promoters. Plays an important role in immune response directly affecting NK maturation and activity, macrophage production of IL12, Th1 development and maturation of CD8+ T-cells. Also implicated in the differentiation and maturation of dendritic cells and in the suppression of regulatory T (Treg) cells development. Its target genes for transcriptional activation activity are: genes invovled in immune response, such as IL7, IL12A/B and IL15, PTGS2/COX2 and CYBB; MHC class I expression, such as TAP1, PSMB9/LMP2, PSME1/PA28A, PSME2/PA28B and B2M and MHC class II expression, such as CIITA. [ref: UniProt]
Pik3cd: PIP3 plays a key role by recruiting PH domain-containing proteins to the membrane, including AKT1 and PDPK1, activating signaling cascades involved in cell growth, survival, proliferation, motility and morphology. Mediates immune responses. Plays a role in B-cell development, proliferation, migration, and function. Required for B-cell receptor (BCR) signaling. Mediates B-cell proliferation response to anti-IgM, anti-CD40 and IL4 stimulation. Promotes cytokine production in response to TLR4 and TLR9. Required for antibody class switch mediated by TLR9. Involved in the antigen presentation function of B-cells. Involved in B-cell chemotaxis in response to CXCL13 and sphingosine 1-phosphate (S1P). Required for proliferation, signaling and cytokine production of naive, effector and memory T-cells. Required for T-cell receptor (TCR) signaling. Mediates TCR signaling events at the immune synapse. Activation by TCR leads to antigen-dependent memory T-cell migration and retention to antigenic tissues. Together with PIK3CG participates in T-cell development. Contributes to T-helper cell expansion and differentiation. Required for T-cell migration mediated by homing receptors SELL/CD62L, CCR7 and S1PR1 and antigen dependent recruitment of T-cells. Together with PIK3CG is involved in natural killer (NK) cell development and migration towards the sites of inflammation. Participates in NK cell receptor activation. Have a role in NK cell maturation and cytokine production. Together with PIK3CG is involved in neutrophil chemotaxis and extravasation. [ref: UniProt]
10 Gene Set Testing
10.1 Molecular Signatures Database
Camera: Performs a competitive gene set test which accounts for inter-gene correlations. Camera evaluates if a set of genes is highly ranked relative to other genes in terms of differential expression.
Barcode plot: This function plots two sets of genes in a ranked list of statistics. Statistics are ranked left to right from smallest to largest. Shaded region in the middle of the plot represents the ranked statistics while the vertical bars reports the positions of the specified subsets. The enrichment worm (line) shows the relative enrichment of the vertical bars in each part of the plot.
Table 5. Significant and selected gene sets from MSigDB libraries
| Gene set | MSigDB | Genes (n) | Direction | P-value |
|---|---|---|---|---|
| LINDSTEDT_DENDRITIC_CELL_MATURATION_B | C2 | 17 | Up | 3.5e-04 |
| KEGG_B_CELL_RECEPTOR_SIGNALING_PATHWAY | C2 | 20 | Up | 3.9e-02 |
| LEE_EARLY_T_LYMPHOCYTE_DN | C2 | 29 | Up | 4.1e-02 |
| SANSOM_APC_TARGETS_DN | C2 | 178 | Up | 4.7e-02 |
| COATES_MACROPHAGE_M1_VS_M2_UP | C2 | 26 | Up | 6.2e-02 |
| COATES_MACROPHAGE_M1_VS_M2_DN | C2 | 22 | Up | 9.6e-01 |
| GSE17186_MEMORY_VS_NAIVE_BCELL_DN | C7 | 37 | Up | 3.2e-02 |
| GSE22886_NAIVE_CD4_TCELL_VS_MEMORY_TCELL_DN | C7 | 31 | Up | 9.8e-03 |
| GSE22886_NAIVE_TCELL_VS_DC_UP | C7 | 53 | Up | 2.1e-02 |
| GSE27786_CD8_TCELL_VS_NKCELL_UP | C7 | 29 | Down | 5.6e-03 |
| GSE29949_MICROGLIA_BRAIN_VS_MONOCYTE_BONE_MARROW_UP | C7 | 88 | Up | 4.7e-02 |
| GSE3982_DC_VS_BCELL_UP | C7 | 55 | Up | 3.5e-02 |
| HALLMARK_INFLAMMATORY_RESPONSE | Hallmark | 64 | Up | 3.4e-02 |
| HALLMARK_INTERFERON_GAMMA_RESPONSE | Hallmark | 63 | Up | 3.8e-02 |
Fig 17. Barcode plot of M1/M2 gene set (up and down) from MSigDB C2 
Fig 17. Figure reports a barcodeplot of gene set COATES_MACROPHAGE_M1_VS_M2_UP (red bars, top of plot) and COATES_MACROPHAGE_M1_VS_M2_DN (blue bars, bottom of plot). The M1/M2 UP is on the border of significance (visible by a one directional enrichment line) while M1/M2 DN is far from significant (visible by non-conclusive enrichment line). Thus there is support for the hypothesis that there is indeed a set of genes related to increase in M1/M2 ratio which are highly ranked in terms of differential expression relative to all other genes.
10.2 Gene Set: Pro-inflammation
Fig 18. Genes specific for pro-inflammation
| GeneSet | Library | NGenes | Direction | PValue |
|---|---|---|---|---|
| Pro-Inflammation | Custom | 7 | Down | 2.6e-05 |
Fig 18. Boxplots reports log-CPM values for genes associated with pro-inflammation which remained post filtering of lowly expressed genes. Genes are a subset of genes form the PCR experiment. Barcodeplot reports the enrichment of genes. Table reports p-value of competitive gene set testing. The boxplots, barcode plot and competitive gen set test confirms that pro-inflammatory genes are downregulated as a group (i.e. in animals which received syngeneic transplantation).
10.3 Gene Set: Microglia
Fig 19. Genes specific for microglia
| GeneSet | Library | NGenes | Direction | PValue |
|---|---|---|---|---|
| Microglia | Custom | 7 | Down | 5.7e-05 |
Fig 19. Boxplots reports genes specific for microglia. Barcodeplot reports the enrichment of genes. Table reports p-value of competitive gene set testing. It is obvious that the expression of microglial genes is lower in mice which received syngeneic IDmBMSCs. This pattern seems to be correlated with the pattern found for pro-inflammatory genes in figure 18. The boxplots, barcode plot and competitive gen set test confirms that pro-inflammatory genes are downregulated as a group (i.e. in animals which received syngeneic transplantation).
11 Comparison with Published Data
Paper: Transcriptome profile of rat genes in injured spinal cord at different stages by RNA-sequencing (Shi et al. 2017)
Authors induced a contusion spinal cord injury in rats and evaluated differential gene expression in animals sacrified at 1, 6 and 28 post injury. Sham animals were included as controls. Differential gene expression of animals at 28 days post injury were used as comparison to data obtained in our study.
Table 6. Comparison of DEG, GO terms and KEGG pathway between Shi et al. 2017 and the current data (injury vs sham)
| Match | DGE.Down | DGE.Up | GO.Term | KEGG.Pathway |
|---|---|---|---|---|
| FALSE | 0.86 | 0.70 | 0.87 | 0.32 |
| TRUE | 0.14 | 0.30 | 0.13 | 0.68 |
| SUM: | 1 | 1 | 1 | 1 |
Table 6. Table reports the fraction of matches (i.e. TRUE is the fraction matched and FALSE the fraction not matched) for differentially expressed genes (up and down), GO terms and KEGG pathway IDs. For the GO terms the number of genes was assumed to be >20 per term (to improve prerequisistes for comparison). Only genes, terms and pathways with p-value <0.05 were considered in the comparison. Considering that different contusion models were used (weight drop vs impactor), difference in time of evaluation (28 days vs 20 days post SCI) and different animal types (rat vs mouse) som significant similarites can be detected.
12 Summary
1. Annotation: Genes are annotated with gene name using their respective ENSEMBL ID.
2. Transformation: Read count matrix is transformed into log-CPM using original library sizes.
3. Filtering: Read count matrix is filtered using log-CPM values (>0 for at least 3 samples).
4. Normalization: Effective library sizes are calculated using the library sizes for the filtered read count matrix and the trimmed mean of M values (TMM) approach.
5. Transformation: Filtered read count matrix is transformed into log-CPM matrix.
6. PCA: Conducted for the 500 genes with highest variance. Proportional variance explained, MDS and loading plots are created.
7. Design matrix: A dummy matrix which indicates which group each sample belongs.
8. Contrast matrix: Contrasts are the group comparisons of interest.
9. Voom transformation: Estimate precision weights for linear modelling to remove dependency between the variance and trhe mean.
10. Linear modelling: Linear modelling using precision weights followed by an empirical Bayes moderation.
11. Differentially expressed genes: Moderated t-statistics are used for determining significantly expressed genes for each contrast. Results are displayed with venn diagrams, mean-difference -and volcano plot and a summary table.
12. Analysis/interpretation: Using hierarchical clustering, heatmap, gene ontology and KEGG enrichment analysis and gene set analysis the difference between the study groups is sought for.
13 Bibliography
[1] R Core Team (2016). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL htts://www.R-project.
[2] Ritchie, M.E., Phipson, B., Wu, D., Hu, Y., Law, C.W., Shi, W., and Smyth, G.K. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research 43(7), e47.
[3] Robinson MD, McCarthy DJ and Smyth GK (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139-140
[4] Law CW, Alhamdoosh M, Su S et al. RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR [version 1; referees: 3 approved]. F1000Research 2016, 5:1408.
14 Setup
This analysis was conducted on:
sessionInfo()## R version 3.3.2 (2016-10-31)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 8.1 x64 (build 9600)
##
## locale:
## [1] LC_COLLATE=Swedish_Sweden.1252 LC_CTYPE=Swedish_Sweden.1252
## [3] LC_MONETARY=Swedish_Sweden.1252 LC_NUMERIC=C
## [5] LC_TIME=Swedish_Sweden.1252
##
## attached base packages:
## [1] grid parallel stats4 stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] pander_0.6.0
## [2] knitr_1.15.1
## [3] pvclust_2.0-0
## [4] boot_1.3-19
## [5] rafalib_1.0.0
## [6] ggrepel_0.6.5
## [7] ellipse_0.3-8
## [8] VennDiagram_1.6.17
## [9] futile.logger_1.4.3
## [10] gplots_3.0.1
## [11] RColorBrewer_1.1-2
## [12] colorspace_1.2-7
## [13] gridExtra_2.2.1
## [14] cowplot_0.7.0
## [15] ggplot2_2.2.1
## [16] data.table_1.10.4
## [17] Mus.musculus_1.3.1
## [18] TxDb.Mmusculus.UCSC.mm10.knownGene_3.2.2
## [19] org.Mm.eg.db_3.3.0
## [20] GO.db_3.3.0
## [21] OrganismDbi_1.14.1
## [22] GenomicFeatures_1.24.5
## [23] GenomicRanges_1.24.3
## [24] GenomeInfoDb_1.8.7
## [25] AnnotationDbi_1.34.4
## [26] IRanges_2.6.1
## [27] S4Vectors_0.10.3
## [28] Biobase_2.32.0
## [29] BiocGenerics_0.18.0
## [30] edgeR_3.14.0
## [31] limma_3.28.21
##
## loaded via a namespace (and not attached):
## [1] bit64_0.9-7 splines_3.3.2
## [3] gtools_3.5.0 assertthat_0.1
## [5] highr_0.6 RBGL_1.48.1
## [7] blob_1.1.0 Rsamtools_1.24.0
## [9] yaml_2.1.13 RSQLite_2.0
## [11] backports_1.0.4 digest_0.6.10
## [13] XVector_0.12.1 htmltools_0.3.5
## [15] plyr_1.8.4 XML_3.98-1.5
## [17] pkgconfig_2.0.1 biomaRt_2.28.0
## [19] zlibbioc_1.18.0 scales_0.4.1
## [21] gdata_2.17.0 BiocParallel_1.6.6
## [23] tibble_1.2 SummarizedExperiment_1.2.3
## [25] lazyeval_0.2.0 magrittr_1.5
## [27] memoise_1.1.0 evaluate_0.10
## [29] graph_1.50.0 BiocInstaller_1.22.3
## [31] tools_3.3.2 stringr_1.1.0
## [33] munsell_0.4.3 lambda.r_1.1.9
## [35] Biostrings_2.40.2 caTools_1.17.1
## [37] RCurl_1.95-4.8 bitops_1.0-6
## [39] labeling_0.3 rmarkdown_1.3
## [41] gtable_0.2.0 DBI_0.5-1
## [43] GenomicAlignments_1.8.4 rtracklayer_1.32.2
## [45] bit_1.1-12 rprojroot_1.1
## [47] futile.options_1.0.0 KernSmooth_2.23-15
## [49] stringi_1.1.2 Rcpp_0.12.7